\chapter{\label{chp:fluo}Time-gated Perturbation Monte Carlo Method for Fluorescence Tomography}
In this chapter, the feasibility of performing FMT with multiple fluorescence compounds based on the time-gated pMC method is investigated. A computationally efficient model that enables the calculation of weight functions for FMT while simulating only the propagation of excitation photons in tissue is derived. This new model is applied to the simultaneous reconstruction of two fluorophores with lifetime contrast in a synthetic mouse model. The information content imparted by single time gates is analyzed in terms of quantification, crosstalk and resolution. Optical reconstructions with multiple time gates are performed to simultaneously image two fluorophores, and the results using different gate selection strategies are compared. Finally, the validity of the model is experimentally established using a slab phantom with objects containing two commercially available fluorophores, and an euthanized mouse with one kind of fluorophore\let\thefootnote\relax\footnotetext{
Portions of this chapter previously appeared as:
\begin{itemize}
\item[] J. Chen, V. Venugopal, and X. Intes, ``Monte Carlo based method for fluorescence tomographic imaging with lifetime multiplexing using time gates,'' \textit{Biomed. Opt. Express}, vol. 2, no. 4, pp. 871-886, 2011.
\item[] V. Venugopal et al., ``Full-field time-resolved fluorescence tomography of small animals,'' \textit{Opt. Lett.}, vol. 35, no. 19, pp. 3189-3191, 2010.
\end{itemize}
}.

\section{Derivation of the Monte Carlo based inverse model for Fluorescence Molecular Tomography}
\noindent In this section, a review of the forward Monte Carlo simulation for fluorescence is provided, based on which the time-resolved inverse model is derived.\\

\noindent {\bf Forward Monte Carlo simulation for fluorescence}
We follow the conventional forward-excitation and forward-emission model outlined by Welch et al~\cite{Welch1997}. The model has been extended to the time domain by first assuming that the emission occurs immediately after the absorption of the excitation photon, then convolving by the lifetime.

In this approach, two Monte Carlo simulations of a complete set of photons are required: one is the propagation of the excitation photons, the other one is the isotropic emission and propagation of the emitted photons. The excitation simulation based on the optical properties at the excitation wavelength $\lambda_x$ is used to calculate the statistical accumulation of the absorbed excitation light $A(r_s, r, t)$, defined as the absorbed photon weights by the fluorophore at position $r$ at time $t$, illuminated by a source at $r_s$. The emission Monte Carlo simulation with the optical properties for the emission wavelength $\lambda_m$ is computed to obtain the emission field $E (r,r_d,t)$ at detector $r_d$ and $t$. Here we define the total absorption coefficients as $\mu_a^x$ and $\mu_a^m$ at $\lambda_x$ and $\lambda_m$, respectively. The total absorption coefficient at $\lambda_x$ is defined as the sum of $\mu_{ab}^x$, the background absorption, and $\mu_{af}^x$, the absorption coefficient contributed by the fluorophore, which is linearly related to the concentration of the fluorophore and the extinction coefficient. Then the effective quantum yield $\eta (r)$ is defined as the probability of a photon to be emitted upon absorption of a photon by the total absorption coefficient:
\begin{equation}
\eta (r) = \frac{\mu_{af}^x (r)\Phi}{\mu_a^x (r)},
\end{equation}
where $\Phi$ is the quantum yield. The time-resolved fluorescence signal without delay caused by the lifetime of the fluorophore is then expressed as:
\begin{equation}
\label{eq:U_f_prime} U'_F (r_s,r_d, t')=\int _\Omega dr^3 W' (r_s,r_d,r,t')\eta
(r),
\end{equation}
where the integration domain $\Omega$ is defined as the entire imaging volume, and the background weight function $W'$ is given by:
\begin{equation}
\label{eq:W_prime} W' (r_s,r_d,r,t')= \int_0^{t'} dt'' A (r_s,r, t'')E
(r,r_d,t'-t'').
\end{equation}
The above non-delay algorithm is summarized in Fig. \ref{fig:fluo_flow}
\begin{figure}[htbp]
\centering
\includegraphics{fluo-1}
\caption{The flowchart of the Monte Carlo based fluorescence simulation.}
\label{fig:fluo_flow}
\end{figure}

Assuming single exponential time-decay for the fluorophore, the fluorescence emission delay can be taken into account by convolving the temporal fluorescence signals with the decay $e^{-t/\tau}$ over time, where $\tau$ is the lifetime of the fluorophore. We can then write the fluorescence intensity measured at $r_d$ and time $t$ for an impulsive excitation at $r_s$ and $t_0=0$ as:
\begin{equation}
\label{eq:lt_conv} U_F (r_s,r_d, t)=\int_0^t dt'U_{F}' (r_s,r_d, t')e^{-
(t-t')/\tau}.
\end{equation}
This general expression is used to compute the fluorescence forward model. In the case of multiple fluorophores, the different fluorescent components can be independently computed and summed to obtain the complex fluorescence time-resolved measurements. In all the \emph{in silico} studies herein, the synthetic fluorescence measurements were based on this approach.\\

\noindent {\bf Calculation of the time-resolved Jacobians}

The forward model described above can be casted into a tomography inverse problem to reconstruct the effective quantum yield $\eta (r)$. In order to compute the time-resolved Jacobians, or the spatial and temporal sensitivity maps with respect to $\eta (r)$, (\ref{eq:lt_conv}) is written as:
\begin{equation}\label{eq:U_f}
U_{F} (r_s,r_d, t)=\int _\Omega dr^3 W (r_s,r_d,r,t)\eta (r),
\end{equation}
where $W$ is the lifetime based weight function defined as:
\begin{equation}
\label{eq:W_conv} W (r_s, r_d, r, t) = \int_0^{t}dt'W' (r_s, r_d, r, t')e^{-
(t-t')\tau},
\end{equation}
Equations~(\ref{eq:W_prime}) and (\ref{eq:W_conv}) give the general expression of the lifetime based weight matrix using Monte Carlo simulations. As (\ref{eq:W_prime}) stands, the background weight function $W'$ can be calculated knowing $A(r_s,r, t)$ and $E(r,r_d,t)$ explicitly. This time convolution is extremely computationally intensive, and simulations with numerous photons at each position $r$ in the region of interest are required to obtain statistically reliable calculation of $E(r,r_d,t)$. Thus this algorithm is computationally inefficient especially in FMT applications where a large number of voxels are considered. However, a more manageable formulation can be derived by considering an assumption commonly employed in FMT, namely that the scattering coefficients are identical at $\lambda_x$ and $\lambda_m$. In the NIR spectral range, the scattering coefficient is expected to vary less than 20\% over the Stokes shift gap of organic fluorophores~\cite{Mourant1998}. Furthermore, under the condition of isotropic scattering, it can be assumed that all fluorescence photons follow the trajectories of the excitation photons after generation~\cite{Liebert2008}. Consider the $i^{th}$ photon in an excitation simulation that propagates from a discrete source point $r_s$ and then detected at the detector position $r_d$ at time $t$. Following the White Monte Carlo (WMC)~\cite{Alerstam2008a} approach, if the photon reaches $r$ at time $t_{s,r}$, the weight of the $i^{th}$ excitation photon $w_i^x$ at $r$ has decreased to:
\begin{equation}
\label{eq:wi_ex} w^x_i (r_s,r,t_{s,r}) = w_{i,0}\exp (-\sum_{j=1}^{p_i}
\mu_{a}^x (r_j) l (r_j)),
\end{equation}
where $w_{i,0}$ is the initial weight of this photon, $r_j(j = 1,...,p_i)$ are the subregions that the photon consequently passes through from $r_s$ to $r$, and $l_i (r_j)$ is the path length that the photon passes at $r_j$. Note that the total number of $p_i$ is photon dependent and the sub-regions for different photons are not necessarily the same. At $r$, the absorbed photon weight is given by:
\begin{equation}
\label{eq:abs_wi} A_i (r_s,r, t_{s,r}) = w_i^x (r_s,r,t_{s,r}) (1-\exp
(-\mu_{a}^x (r)l_i (r))).
\end{equation}
The absorbed weight multiplied by $\eta (r)$ is then converted to the initial weight of the fluorescence photon. We can calculate the final weight of the fluorescence photon detected at $r_d$ as:
\begin{equation}
\label{eq:final_w} w_i^m (r,r_d, t_{r,d}) = A_i (r_s,r, t_{s,r})\eta (r)\exp
(-\sum_{j=p_i+1}^{q_i}\mu_{a}^m (r_j)l_i (r_j)),
\end{equation}
where $t_{r,d}$ is the time that the photon passes from $r$ to $r_d$ and $r_j(j = p_i+1,...,q_i)$ is the photon path from $r$ to $r_d$. Note that sum of $t_{s,r}$ and $t_{r,d}$ is the total time $t$ that the excitation photon travels from $r_s$ till is detected, and that the subregions $r_j(j = 1,...,q_i)$ denote the total photon path. By summing up $n$ detected emission photons for $r_s$ and $r_d$, we can obtain $U'_F$ defined in (\ref{eq:U_f_prime}). Comparing (\ref{eq:final_w}) to (\ref{eq:U_f_prime}), we have:
\begin{equation}
\label{eq:W_prime_m} W' (r_s, r_d, r, t)= \sum_{i=1}^n A_i (r_s,r, t) \exp
(-\sum_{j=p_i+1}^{q_i}\mu_{a}^m (r_j)l_i (r_j)).
\end{equation}
Inserting (\ref{eq:wi_ex}) and (\ref{eq:abs_wi}) to (\ref{eq:W_prime_m}), then we have:
\begin{equation}
\begin{split}
\label{eq:W_prime_1} W' (r_s, r_d, r, t)= \sum_{i=1}^n w_{i,0}\exp
(-\sum_{j=1}^{p_i}\mu_{a}^x (r_j)l_i (r_j))
\times (1-\exp (-\mu_{a}^x (r)l_i (r))\\
\times \exp (-\sum_{j=p_i+1}^{q_i}\mu_{a}^m (r_j)l_i (r_j)).
\end{split}
\end{equation}
According to (\ref{eq:W_prime_1}), the background weight matrix for excitation source $r_s$ and detector $r_d$ at time $t$ can be easily calculated using the photon paths and the absorption coefficients. This formulation should be employed in the case that the absorption coefficient at the excitation significantly differs from that at the emission wavelength, which might be the case for wavelengths below 700nm. However, the absorption optical spectra of bio-tissues is relatively flat in tissue in the NIR window~\cite{Roggan1995}. Thus for simplicity, and following standard assumption in FMT, we assume that the absorption coefficients at the excitation and emission wavelengths are identical, that is, $\mu_a^x = \mu_a^m$. We therefore obtain:
\begin{equation}
\label{eq:W_prime_2} W' (r_s, r_d, r, t)= \sum_{i=1}^n w_i^x (r_s,r_d,t)
(1-\exp (-\mu_{a}^x (r)l_i (r)),
\end{equation}
where $w_i^x (r_s,r_d,t)$ is the final weight of the $i^{th}$ excitation photon from $r_s$ and detected by $r_d$ at $t$, defined as:
\begin{equation}
\label{eq:w_final} w_i^x (r_s, r_d, t)= w_{i,0}\exp (-\sum_{j=1}^{q_i}\mu_{a}^x
(r_j)l_i (r_j)).
\end{equation}
If $\mu_{af}^x (r)l_i (r)$ is small, which is the case in the voxelized geometry, from Taylor expansion, we can further simplify (\ref{eq:W_prime_2}) to:
\begin{equation}
\label{eq:W_prime_3-1} W' (r_s, r_d, r, t)= \sum_{i=1}^n w_i^x
(r_s,r_d,t)\mu_{a}^x (r)l_i (r).
\end{equation}
Integrating (\ref{eq:W_prime_3-1}) over time results in a CW reconstruction algorithm similar to the method proposed by Zhang et al~\cite{Zhang2009}, that the photon paths weighted by their final fluence rate seen by the detector is the sensitivity function. Note also that the background weight matrix for fluorescence calculated by (\ref{eq:W_prime_3}) is precisely the weight function for absorption perturbations. This indicates that the absorption and fluorophore distribution can be estimated by using the same weight matrix calculated by one excitation Monte Carlo simulation. Moreover, it is noteworthy that if there are different fluorophore species existing in the region of interest, only one excitation simulation is required to obtain the weight matrices $W_k$ ($k=1,...,N_F$, where $N_F$ is the number of species) for different lifetimes $\tau_k (k=1,...,N_F)$ according to (\ref{eq:W_conv}). This allows for a fast and efficient computational implementation to simultaneously reconstruct multiple fluorophores based on one forward simulation. Note that for simplicity, all the above equations use $r_s$ and $r_d$ as discrete points, but they can be extended to any illumination or detection strategies applied to complex boundaries by collecting the photons generated or detected inside the areas. In this study, we employed the widefield illumination strategy discussed in the previous Chapter.\\

\noindent {\bf Inverse model for fluorescence tomography}

To formulate the inverse problem, we employed a Born formulation in which the emission field is normalized by the excitation field~\cite{Soubret2005} by dividing the experimental time-domain emission measurements with the CW excitation flux at the same position. We therefore have:
\begin{equation}
\label{eq:normalize} \frac{M^m (r_s,r_d, t)}{M^x (r_s,r_d)}=\frac{ \alpha}{U^x
(r_s,r_d)}\sum_{k=1}^{N_F}\int d^3rW_k (r_s,r_d, r,t)\eta_k (r),
\end{equation}
where $\alpha$ incorporates unknown constants associated with wavelength-dependent gains and attenuations that can be measured once for every imaging system, $M^m (r_s, r_d, t)$ is the total signal from all fluorophores with different lifetimes at the emission wavelength measured at the boundary position $r_d$ and $t$ excited from the source at $r_s$, $M^x (r_s, r_d)$ and $U^x (r_s, r_d)$ are the measured and simulated, respectively, total excitation flux measured by a photon-detector at $r_d$. This normalization efficiently mitigates the dependence of the detected fluorescent signal on the optical properties of the examined tissue~~\cite{Soubret2005, Vinegoni2009}, thus we did not model the absorptive heterogeneities associated with the different organs in our synthetic murine model.

The normalized fluorescence signals corresponding to different source-detector pairs and time gates (left term in (\ref{eq:normalize})) can be stacked up as a vector $\gamma\in R^{N_m}$ ($N_m$ is the total number of measurements, that is, the product of the number of source-detector pairs and the number of time gates) where it can be used as part of an inverse problem. Similarly, the $k^{th}$ effective quantum yield for the entire imaging volume can be written as $\eta_k^\Omega = [\eta_k(r_1),...,\eta_k(r_{N_V})]^T$, where $N_V$ is the total number of voxels in the region of interest. Overall, the forward problem can be expressed in a matrix form:
\begin{equation}
\left[\begin{array}{cccc}
\gamma_1 \\
\vdots \\
\gamma_{N_m}
\end{array}\right]
 =
\left[\begin{array}{ccc}
\beta_1W_{1,1}^\Omega & \dots & \beta_1W_{1,N_F} ^\Omega \\
\vdots & \ddots&\vdots \\
\beta_{N_m}W_{N_m,1}^\Omega & \dots & \beta_{N_m}W_{N_m,N_F}^\Omega
\end{array}\right]
\left[\begin{array}{c}
\eta_1 ^\Omega \\
\vdots \\
\eta_{N_F}^\Omega
\end{array}\right],
\end{equation}
where $W_{s,k}^\Omega$ is the weight function for the $s^{th}(s = 1,...,N_m)$ measurement and the $k^{th}$ fluorophore species over the volume $\Omega$, and $\beta_s = {\alpha_s}/{U^x_s}$ is the normalization factor for the $s^{th}$ measurement and the corresponding weight function. The above equation has a general form $Ax=b$, then $x = [\eta_1^\Omega,...,\eta_{N_F}^\Omega]^T$ can be computed by solving this linear system of equations. Notice here the reconstructions for all fluorescent species provided herein are simultaneous based on monochromatic data and the measurements employed are direct measurements that are not derived from a pre-processed unmixing algorithm. Because of the ill-posedness of the inverse problem, the reconstruction method in fluorescence tomography usually sets up a least square optimization problem minimizing an objective function given by:
\begin{equation}
\label{eq:obj_func} \Psi = \| Ax-b\|^2+\|\lambda (r)x\|^2,
\end{equation}
where $\lambda (r)$ is a spatially variant regularization parameter to compensate for the spatial dependence of the contrast and resolution in the reconstruction~\cite{Pogue1999}.

%flowchart and some calculation details
%\section{\label{sec:flow}Reconstruction algorithm}
In summary, our approach to reconstruct the effective quantum yield is schematically depicted in Fig.~\ref{fig:fluo-flow}. The capital letters shown in the square brackets below refer to the individual steps in this scheme. [A] represents a time-resolved Monte Carlo simulation for excitation photons, following the WMC approach to analytically compute the light transport in tissue. This simulation can be adapted to different illumination and detection strategies. To calculate the background weight matrix in [B], we directly add the weighted paths to the weight matrix along with the Monte Carlo simulation, instead of storing all the photon trajectories. This reduces the post processing time for weight matrix calculation as well as the storage space for the huge number of photon paths. Note that any of (\ref{eq:W_prime_1}), (\ref{eq:W_prime_2}) and (\ref{eq:W_prime_3-1}) can be used in [B]. It should be noted that (\ref{eq:W_prime_3-1}), wherein the absorption coefficient at $\lambda_x$ and $\lambda_m$ are assumed to be identical, can be used for computational efficiency. The lifetime-based weight matrices are then calculated by (\ref{eq:W_conv}) in [C]. Knowing the experimental measurements, a conjugate gradient method is applied in [D] to minimize the objective function defined in (\ref{eq:obj_func}) and obtain the fluorescence yield distribution. The fluorescence measurements in this chapter are either from an actual experiment or an emission Monte Carlo simulation at $\lambda_m$ based on the absorption probability calculated in [A].
\begin{figure}[htbp]
\centering
\includegraphics{fig_flow}% Here is how to import EPS art
\caption{\label{fig:fluo-flow} Block diagram summarizing the steps to perform fluorescence reconstruction.}
\end{figure}

\section{\textit{In silico} fluorophore lifetime multiplexing study}
\subsection{\label{sec:sim_setup}Simulation setup}
We first validated the proposed algorithm \emph{in silico} in a model replicating small animal fluorescence imaging. The Monte Carlo simulations were performed on a mouse geometry with the kidneys and the skin shape extracted from a whole-body atlas~\cite{Dogdas2007} (Fig.~\ref{fig:sim_setup}a). The entire field of view consisted of $ 89 \times 33 \times 19$ voxels with size $1$mm$^3$. The optical properties were set to $\mu_{ab}^x = 0.3$cm$^{-1}$, $\mu_s'^x= 25$cm$^{-1}$, $\mu_{a}^m = 0.36$cm$^{-1}$, $\mu_s'^m= 20$cm$^{-1}$ over the entire body, to simulate the different optical properties of mouse tissues at different wavelengths in the NIR window \cite{Roggan1995}. The kidneys were considered to be labeled with two distinct fluorophores with lifetime of 0.5ns (black) and 1ns (red), respectively. The $\mu_{af}$ for both kidneys were set to 0.1 times of the background absorption coefficient $\mu_{ab}^x $, and the quantum yield was set to 1. 36 widefield illumination patterns (Fig.~\ref{fig:sim_setup}b) and a grid of 64 point detectors in transmission geometry were arranged covering a 31mm $\times$ 27mm $\times$ 18mm volume spanning the abdomen. The sources and detectors are projected to the surface along the $z$ axis. The initial positions of the photons were randomly spanned (that is, uniformly distributed) over the illuminated area to accommodate the broad field illumination. The detectors had a 2mm separation along the $y$ axis and a 2.2mm separation along the $x$ axis with a radius of 1mm.
\begin{figure}[htbp]
\centering
\includegraphics[width=4.5in]{fluo-4}% Here is how to import EPS art
\caption{\label{fig:sim_setup}(a) The mouse model geometry and arrangement of sources (top) and detectors (bottom) used for the simulations. The fluorescent kidneys are marked as red and black, according to different lifetimes. (b) The illumination patterns (red: on, black: off).}
\end{figure}

In the Monte Carlo simulation, $10^9$ photons were consecutively launched for each pattern source on a massively parallel environment (blue gene, CCNI at RPI) to calculate the excitation and fluorescence measurements as described in Section~\ref{sec:meth_tr_mc}. Under the assumption that the optical properties are equal at $\lambda_x$ and $\lambda_m$ in (\ref{eq:W_prime_3}), the time-resolved jacobians for FMT were simultaneously computed by convolving the excitation field jacobians by the lifetime decay. The temporal profile was recorded over a 6ns time window with 20ps resolution and a gate-width of 200ps (total of 120 gates). These parameters were selected to replicate the operating characteristics of our preclinical imaging platform described in Chapter~\ref{chp:pmc} and correspond to the typical minimum gate width. Note that we do not explicitly include noise in the simulation data. However, the Poisson noise inherent in time-resolved studies also exists in Monte Carlo simulations. As the excitation measurements and fluorescence measurements are both directly simulated in this work (not computed through the multiplication of a Jacobian with a synthetic phantom), both measurements have an uncorrelated Poisson noise that replicates experimental conditions.

The average calculation time for an excitation or emission MC simulation for one pattern source and all detectors and gates was less than 8 minutes (on 1024 nodes). Fig.~\ref{fig:raw_data} shows typical simulated fluorescence signals and the background weight matrices calculated using (\ref{eq:W_prime_3-1}) at different time gates.
\begin{figure}[htbp]
\centering
\includegraphics[width=4in]{fluo-5}
\caption{\label{fig:raw_data} (a) Representative simulated signals from the fluorophores with $\tau_1 = 0.5$ns~(blue) and $\tau_2 = 1$ns~(red), and the total fluorescence signal (black) used in the reconstructions. All the signals are generated using the mouse phantom and a single transmittance S-D pair as shown in (b) and (c). (b) Weight matrices for $\tau = 1$ns at the 3 gates (gray) presented in (a). (c) The direct detector readings from the region of interest (black boxes) from the simulation on the mouse phantom. The bars in (b) and white boxes in (c) indicate the illumination pattern and the crosses in (b) and (c) indicate the detector used in (a).}
\end{figure}

\subsection{\label{sec:sim_single}Gate information content analysis}
In order to quantitatively assess the information content at single gates, we need to establish quantitative performance metrics. We investigated the quantification, crosstalk and resolution based on the reconstructions using one gate from gate 10 to gate 145 at a separation of 5 gates. The quantification of the reconstructed $\eta_1$ is defined as:
\begin{equation}
Q_1=\frac{\max[\eta_1(\Omega_1)]}{H_1(\Omega_1)}\times 100\%,
\end{equation}
where $\Omega_1$ denotes the known location of the inclusions with lifetime $\tau_1 = 0.5$ns and $H_1$ denotes the expected value of $\eta_1$. The crosstalk is defined as:
\begin{equation}
X_1 =
\frac{\max[\eta_2(\Omega_1)]}{\max[\eta_1(\Omega_1)]\times 100\%},
\end{equation}
to quantify the separability of the two inclusions with lifetime contrast. The quantification and yield crosstalk for the 1ns lifetime component was similarly evaluated. The spatial resolution is measured as the full width at half maximum (FWHM) of the point spread functions (PSF). PSFs are created by simulating a perturbation at a single point. The simulated perturbation, $x_{sim}$, is a vector of zeros except for one target pixel at the center of the reconstructed volume with a value of one. We then generated simulated measurements $y_{sim}=Ax_{sim}$ and reconstructed the image. The cube root of the volume of voxels within half the peak reconstructed value is defined as FWHM. Plots of the quantification, crosstalk and FWHM of PSFs are shown in Fig. ~\ref{fig:gate_analysis}. Note that for all these reconstructions, the inverse problem size was identical and that the same iteration number (150) was used in the CG algorithm. This ensured consistency between the different reconstructions.
\begin{figure}[htbp]
\centering
\includegraphics[width=\textwidth]{fluo-6}
\caption{\label{fig:gate_analysis}Information content evaluation based on single gate reconstructions in terms of (a) quantification, (b) crosstalk and (c) resolution.}
\end{figure}

From the plots of performance metrics, it is clear that for the reconstructed object with the shorter lifetime, the quantitative accuracy is decreasing with time while the crosstalk is increasing. Conversely, for the reconstructed object with the longer lifetime, an opposite behavior is observed with equivalent quantification and crosstalk around the maximum gate. It is noteworthy that the reconstruction using the latest gate ($t = 2.8$ns) results in the most accurate quantification and the least crosstalk for the compound with the longer lifetime, implying the necessity of the late gates for effective separation of the objects with lifetime contrast. It is impossible to achieve the best quantification and minimum crosstalk using only the early gates because of the comparable contribution of the two fluorophores at the early time points (Fig.~\ref{fig:raw_data}). As expected, the steadily increasing FWHM of PSFs with time indicates that the early gates provide improved resolution when compared to the late gates \cite{Soloviev2009, Niedre2010}. The better resolution at early gates is due to the narrower probing volume by photons detected at the early time points, as seen in Fig.~\ref{fig:raw_data}b. However, there is no significant reduction in the resolution of the reconstructed yield distributions for the two lifetimes when using gates later than the maximum gate. The reconstruction at the 10$^{th}$ gate ($t=0.2$ns) does not follow the trend due to the poor statistics of the weight matrix at the early gates, where the number of detected photons is insufficient to generate reliable statistics. It results in increased artifacts in the reconstructions when using very early gates. From the above analysis, we conclude that multiple fluorophores with lifetime contrast cannot be efficiently separated with minimal crosstalk if a single gate of the TPSF is used.


\subsection{\label{sec:sim_rec}Lifetime multiplexing tomography based on multiple gates}
\noindent To utilize the various information carried by different gates, we investigated reconstructions using multiple gates simultaneously. In this regard, 5 sets of time-gated measurements consisting of 5 gates each were considered. The number of gates was limited to 5 to reduce the memory burden as all the gates are simultaneously used in one inverse problem. The 5 different sets simulated in this study are illustrated in Fig.~\ref{fig:comb_recon}a. Sets 1, 3, 5 have gates spaced at uniform intervals covering different time ranges, corresponding to the early to the maximum gate, full span, and the maximum to the end of the TPSFs. Sets 2 and 4 span the full time range, with non-uniform spacings exponentially scaled towards the rising and the decaying edges of the TPSFs respectively. This is a heuristic but straightforward approach to investigate the reconstruction performance using multiple combined gates.
\begin{figure}[htbp]
\centering
\includegraphics{fluo-7}
\caption{\label{fig:comb_recon} (a) Sets of time gates investigated in multi-gate reconstructions. The reconstructions using combined gates of set 1 and set 5 are displayed in (b) and (c), respectively.}
\end{figure}

Table~\ref{tab:table1} summarizes the quantitative metrics described in the previous section for the 5 sets investigated. It should be noted that the introduction of late gates improves the quantification of both fluorescence yields. Moreover, significantly less crosstalk than any of the single-gate reconstructions was achieved by using combined gates. The reconstruction using early gates (Set 1) has superior resolution compared to the reconstruction using late gates. However, the crosstalk for this case, is significantly higher than those applying late gates. For sets 2-4, the quantification and crosstalk are improving with the gate selection shifting to late gates, with similar performance in resolution due to the use of an early gate. It should be noted that the most accurate quantification and minimum crosstalk is obtained for both lifetimes when using the maximum and late gates (Set 5). However, without applying any early gates in reconstruction, the resolution worsened by 20\% compared to the other cases.

In order to ascertain the effect of assuming equal optical properties at excitation and emission wavelengths, a second set of time-resolved Jacobians were computed for the same set of gates using (\ref{eq:W_prime_1}) (A-E). The reconstruction results when using sets are also given in table~\ref{tab:table1}. Comparing sets 1-5 and sets A-E, it is determined that the assumption of equal optical properties at excitation and emission wavelengths results in less than 5\% quantification error for all three criteria.

\begin{table}
\caption{\label{tab:table1}{Quantitative comparison of the multi-gate reconstructions.}}
\begin{center}
\def\temptablewidth{0.7\textwidth}
{\rule{\temptablewidth}{1pt}}
\begin{tabular}{ccccccc}
Set & \multicolumn{2}{c}{Quantif. (\%)}&\multicolumn{2}{c}{Crosstalk (\%)}& \multicolumn{2}{c}{FWHM (mm)} \\
&$\tau_1$&$\tau_2$&$\tau_1$&$\tau_2$&$\tau_1$&$\tau_2$\\
\hline
1 & 44.53& 45.65& 39.74& 50.43 & 2.33& 2.44\\
2 & 45.89& 49.22& 35.69& 36.50 & 2.52& 2.64\\
3 & 46.98& 52.19& 32.25& 33.42 & 2.69& 2.78\\
4 & 49.04& 53.96& 29.20& 30.66 & 2.74& 2.85\\
5 & 49.77& 63.05& 24.92& 24.10 & 3.21& 3.16\\
\hline
A & 47.64& 45.72& 37.19& 49.19 & 2.19& 2.32\\
B & 48.59& 50.21& 33.34& 33.98 & 2.43& 2.46\\
C & 49.35& 53.40& 31.01& 31.45 & 2.63& 2.67\\
D & 50.47& 54.39& 27.34& 30.05 & 2.72& 2.76\\
E & 55.76& 66.45& 21.77& 23.19 & 3.17& 3.11
\end{tabular}
{\rule{\temptablewidth}{1pt}}
\end{center}
\end{table}

%\begin{table}
%\caption{\label{tab:table2}Quantitative comparison of the multi-gate
%reconstructions \emph{corrected mua using eq. 11... slightly better}}
%\centering
%\begin{tabular}{ccccccc}
%Set & \multicolumn{2}{c}{Quantif. (\%)}&\multicolumn{2}{c}{Crosstalk (\%)}& \multicolumn{2}{c}{FWHM (mm)} \\
%&$\tau_1$&$\tau_2$&$\tau_1$&$\tau_2$&$\tau_1$&$\tau_2$\\
%\hline
%1 & 47.64& 45.72& 37.19& 49.19 & 2.19& 2.32\\
%2 & 48.59& 50.21& 33.34& 33.98 & 2.43& 2.46\\
%3 & 49.35& 53.40& 31.01& 31.45 & 2.63& 2.67\\
%4 & 50.47& 54.39& 27.34& 30.05 & 2.72& 2.76\\
%5 & 55.76& 66.45& 21.77& 23.19 & 3.17& 3.11
%\end{tabular}
%\end{table}

\section{Experimental study}
\subsection{Phantom experiment}
\noindent
The lifetime based Monte Carlo algorithm was experimentally validated by using the time-domain small animal molecular imaging system described in Chapter~\ref{chp:pmc}. 36 binary sliding-bar patterns are reflected on the programmable DMD chip and projected on the imaging chamber to provide an illumination area of 35mm $\times$ 20mm (matching the dimensions of a small animal torso). The transmitted signal is detected with a gate-width of 300ps and a time resolution of 40ps. Photon measurements at the excitation and fluorescence wavelengths were acquired consecutively using bandpass filters. The incident fields for excitation and the fluorescence data for excitation were collected at 780nm and 832nm, respectively. The laser power was independently optimized for both the excitation and fluorescence measurements to maintain a maximum signal level of 4000 counts. Further, the acquired images are post-processed to generate 1mm $\times$ 1mm detector measurements by averaging over 7$\times$7 pixels.

%\subsection{\label{sec:exp_phantom}Phantom setup}

A polycarbonate slab phantom (Fig.~\ref{fig:expe}b) that consisted of two tubes having an inner diameter of 1.5mm was used to experimentally assess the performance of the algorithm. To achieve lifetime contrast, tube I, II were filled with 2.5nmol of Indocyanine Green ($\tau_1 = 0.45$ns) and IRDye 800CW ($\tau_2 = 0.8$ns) dissolved in 180$\mu$L ethanol, respectively. The effective quantum yields of the two tubes was externally calibrated and tube I was found to have 1.5 times the yield of tube II. The tank was filled with a mixture of Intralipid-20\% (Sigma-Aldrich) and a water-soluble NIR Dye (Epolight 2717, Epolin) diluted in water to simulate the average optical properties ($\mu_a = 0.16$cm$^{-1}$ and $\mu_s' = 17$cm$^{-1}$) of murine models.
\begin{figure}[htbp]
\centering
\includegraphics[width=2.5in]{fluo_phantom_setting}
\caption{\label{fig:expe} Phantom setup. The red box shows the pattern illumination area. All dimensions are in mm.}
\end{figure}

%\subsection{\label{sec:exp_rec}Experimental reconstruction}

For the MC computation, the phantom was digitized into a 28 $\times$ 40 $\times$ 14 volume image with 1mm $\times$ 1mm $\times$ 1mm voxels. A background weight matrix as defined in (\ref{eq:W_prime_3-1}) covering the field of view of 24 $\times$ 34 $\times$ 14 area, resulting in a total of 11,424 voxels, was calculated. The accuracy of the inverse model with experimental data is affected by the heterogeneous spatial distribution of the experimental pattern on the phantom. The illumination spatial heterogeneity is taken into account in the Monte Carlo simulation by scaling the initial photon weight to the illumination intensity at the injection position using experimental calibration patterns. The system impulse response functions (IRF) for all sources was convolved into the model before tomographic inversion, to match the model prediction with the experimental measurements. Since the lifetimes are known in advance in this controlled phantom experiment, the weight functions are generated using the shorter lifetime at $\tau_1 = 0.45$ns (ICG) and the longer lifetime at $\tau_2=0.8$ns (IRdye\texttrademark) using (\ref{eq:W_conv}). In \emph{in vivo} applications, the lifetime of the dye within a control animal should be separately estimated due to the possible fluctuation of the fluorophore lifetime under different micro environments.

To achieve better separation between the two fluorophores while maintaining resolution, we selected the measurements at 50\% of the rising edge, the maximum, 75\%, 50\%, 25\% of the decaying edge of the convolved TPSFs (mimicking gate set 4) to recover the fluorescence yield distribution based on the analysis of the simulated result. Gates with less than 400 photons were discarded in the reconstruction process to maintain an appropriate signal- to-noise ratio (SNR). For comparison, we also reconstructed the fluorophore distribution using CW data. The CW weight functions for the inversion were calculated by integrating the corresponding TD weight functions over time and the TPSFs were integrated over time to obtain the CW datatype.
\begin{figure}[htbp]
\centering
\includegraphics{fluo-8}
\caption{\label{fig:fluo-exp_result} Phantom experiment results using (a) TD data and (b) CW data. The images are normalized to the maximum of the reconstruction.}
\end{figure}

The reconstruction results are shown in Fig.~\ref{fig:fluo-exp_result}, similar as the 3D visualization of the \emph{in silico} results. The two objects are separated with crosstalk $X_1 = 31.23\%$ and $X_2 = 28.77\%$. The ratio of the mean reconstructed yields within the 50\% isovolume of tube I and II was found to be 1.53. The average dimension of the reconstructed tubes was 37.5\% larger along the $x$~axis and 87.5\% larger along the $z$~axis~(depth). This difference in resolution is due to the transmittance arrangement of the source-detector pairs, which smears the resolution along the $z$~axis. Hence, the diameter of the reconstructed isosurface is approximately 2 times bigger than the actual diameter due to the loss in the depth dimension, which can be attributed to the high scattering, and to the position of the two objects in the middle of the phantom, where the weight matrix has the least sensitivity. The CW reconstruction presented in Fig.~\ref{fig:exp_result}b, shows poor resolution (103\% increase in depth) and no separation (100\% crosstalk).

\subsection{\label{sec:fluo_mouse_exp}Animal experiment}
\noindent In a second experiment, we investigated the applicability of the algorithm in a small animal with complex boundaries. The experiment is performed on a freshly euthanized mouse with a 13mm long with 1mm diameter placed in the thoracic cavity. The tube contains 14$pmol$ of IRDye 800CW R(LiCOR) in 1$\mu L$ ethanol. The animal was restrained in the imaging chamber and consecutively imaged on a $\mu$CT scanner (Scanco Viva CT40) and the optical tomography platform. The CT scans were used to extract the surface geometry and pinpoint the localization of the object (Fig. \ref{fig:fluo_mouse_setting}a).

The tomographic imaging protocol used 64 bar patterns projected over a 33mm $\times$ 18mm area across the upper torso of the animal (Fig. \ref{fig:fluo_mouse_setting}b and \ref{fig:fluo_mouse_setting}c). 97 point detectors measured the excitation and fluorescence photons transmitted through the animal. The mouse model extracted from the CT images was discretized into 0.5mm$\times$0.5mm$\times$0.5mm voxels. The optical measurements were registered with this model using fiducial markers on the imaging chamber. The 25\% rising TG is selected as the datatype and the corresponding weight functions are computed using the average background properties of the small animal thoracic cavity ($\mu_a$ = 0.3cm$^{-1}$ and $\mu_s'$ = 25cm$^{-1}$) \cite{Niedre2010}.
\begin{figure}[htbp]
\centering
\includegraphics{fluo_mouse_setting}
\caption{\label{fig:fluo_mouse_setting} (a) The segmented CT image and the position of the tube in the chest cavity; the boundary extracted from CT-scanned images overlayed with (b) half-space illumination patterns and (c) normalized Born measurements. The red box and the black circles show the entire illuminated area and the detector positions, respectively.}
\end{figure}

Figs. \ref{fig:fluo_mouse_result}b-c show the transverse and coronal cross-sections of the reconstructed volume. The reconstructed inclusion had a length of 10.8mm and 3mm diameter. In this experiment the fluorescent inclusion is accurately localized in the chest cavity of the animal cadaver demonstrating the feasibility to apply the Monte Carlo method for FMT in small animals. The tomographic measurements of the excitation and fluorescence photons were completed in less than 30 minutes using the widefield strategy. These strategies can be easily extended to whole-body multispectral FMT allowing efficient separation and 3D reconstruction of multiplexed exogenous fluorophores based on a spectral and temporal contrast.
\begin{figure}[htbp]
\centering
\includegraphics{fluo_mouse_result}
\caption{\label{fig:fluo_mouse_result}(a) Reconstructed 3D Volume with the CT scan. (b) Coronal slice of the reconstructed volume at z = 6.5mm. (c) Transverse slice of the volume at $y =$21.5mm.}
\end{figure}

%\section{Optimum gate selection for robust multi-lifetime tomography}
%\noindent
%In tomography settings, if a large number of source-detector pairs or a region of interest including many voxels is employed, using more gates is leads to a larger memory taken and thus not practical. Moreover, if the optimum gate set can be predicted before measuring the fluorescence, it will take much shorter time for the fluorescence detection, which is the most time consuming part of the experiment. Based on these consideration, it is requisite to have the least number of gates that can provide similar level of accuracy in reconstructions.
%
%In this section, an approach based on error propagation is developed to choose best gates combination for simultaneous reconstruction of fluorephores with different lifetime. This method is applied to the mouse model with two kidneys and compared with the previous results with empirically selected gates.

%\subsection{Solution for double lifetime multiplexing}
%Consider two lifetimes $\tau_1$ and $\tau_2$, for any 2 gates with measured intensity $I_1$ and $I_2$ we have
%\begin{equation}\label{eq:optm-3}
%\begin{split}
%  I_1 = k_1A_1 + k_2 B_1,\\
%  I_2 = k_1A_2 + k_2 B_2.
%\end{split}
%\end{equation}
%
%where $A_n$ and $B_n$ are the exponential decay component for $\tau_1$ and $\tau_2$ and $k_n$ is the coefficients for them. The solution of the component coefficients can be then calculated as
%\begin{equation}\label{eq:optm-2}
%\begin{split}
%k_1 = \frac{B_2 I_1 - B_1 I_2} {A_1 B_2 - A_2 B1},\\
%k_2 = \frac{A_1 I_2 - A_2 I_1}{A_1 B_2 - A_2 B_1}.
%\end{split}
%\end{equation}
%
%When attempting to use more time points than the lifetimes, the method of least square, as a standard approach for overdetermined systems, can be employed to obtain the approximation solution. Assuming there are $N$ gates available, we have
%\begin{equation}\label{eq:optm-3}
%\begin{split}
%  I_1 = k_1A_1 + k_2 B_1,\\
%  I_2 = k_1A_2 + k_2 B_2,\\
%  \vdots\\
%  I_N = k_1A_N + k_2 B_N.
%\end{split}
%\end{equation}
%
%If $n$ gates are selected from the total $N$ gates for estimating $k_n$, from the least square solution for this linear system, we have
%\begin{equation}\label{eq:optm-1}
%\begin{split}
%y_i = \frac{I_i}{B_i},\\
%t_i = \frac{A_i}{B_i},\\
% \bar{y} =\frac{1}{n}\sum y_i,\\
% \bar{t} = \frac{1}{n}\sum t_i,\\
%  k_1 = \frac{\sum t_i y_i -n \bar{t} \bar{y}}{\sum t_i^2 - n \bar{t}^2},\\
% k_2 = \bar{y}-k1 \bar(t).
%\end{split}
%\end{equation}
%% C = \sum{t_i^2}-n\bar{t}^2,\\
%% D = \sum t_i y_i-n \bar{t} \bar{y},\\
%where i is the $i^{th}$ gate out of $N$.
%
%
%\subsection{Error propagation and optimum gates for tomography}
%Propagation of error describes the effect of variables' error on the error of a function based on them. The uncertainties of the desired quantity can be calculated using the following formula based on the changing variables (usually experimental measurements) that have uncertainties due to measurement limitations such as instrument precision.
%\begin{equation}\label{a}
%  \sigma = \sqrt{\left(\frac{\partial f}{\partial x_1} \right)^2 \sigma_1^2+\left(\frac{\partial f}{\partial x_2} \right)^2 \sigma_2^2+\dots+ \left(\frac{\partial f}{\partial x_n} \right)^2 \sigma_n^2}
%\end{equation}
%
%This method is usually used when the value of the quantity of interest is not easy to measure directly but can be expressed by other quantities that are easy to measure.
%
%In FLIM applications, the method of error propagation is mostly applied to determine optimal gating schemes for accurate lifetime estimation. Since optical measurements have a Poisson type of noise, that is, the standard deviation of the measurements is equal to the square root of the expected value, the error of the measurements can readily propagates to the error of the estimation of the lifetimes, once the lifetime calculation is a function of the measurements.
%
%Now we apply this method for optimizing the gates selection for time-gated fluorescence molecular imaging (or FRET). For simplicity, the lifetime is treated as a priori, which can be found accurate in literature, or from a global fit of the fluorescence images. Due to this assumption, the time-gated profile at every detector for each lifetime can be formed separately. The light intensity measurement, at each time point, becomes a linear combination of these profiles. If no noise is considered, for a given number of lifetimes, the unknowns, which are the coefficient for each component, can be solved directly using the same number of measurement points. It is always more accurate when more time points (gates) are used for fitting.
%
%Based on the error propagation formula, the propagated error on k1/k2 can be calculated for any gate combination. An explicit expression for the propagated error for $k1$ and $k2$ was derived and coded to avoid the slow symbolic calculation in Matlab. To directly calculate error is 200 times faster than to use the public code based on symbolic written. The results using a publicly available code for calculating propagated error is compared with this code to verify the formula.
%
%
%%The propagated error of $k_1/k_2$ is evaluated to access the accuracy of both $k_1$ and $k_2$.
%
%%The detailed calculation of the error of the k1/k2 ratio is as following
%%\begin{equation}\label{eq:optm-1}
%%\begin{split}
%% \bar{y} =\sum(I./B1(i,:))/n,\\
%% \bar{t} = \sum(A1(i,:)/B1(i,:))/n,\\
%% C = \sum((A1(i,:)./B1(i,:))^2)-n*t_bar.^2,\\
%%  D = \sum(A1(i,:).\times I./B1(i,:).^2)-n\times t_bar \times y_bar,\\
%%  ratio_std(i) = sum((C/D^2./B1(i,:).*(D/n-y_bar*(A1(i,:)./B1(i,:)-t_bar))).^2.*I);
%%\end{split}
%%\end{equation}
%
%\subsection{\emph{In silico} results}
%
%We first applied the formula on pure decay curves for comparing with the result from the FLIM theory. Here the shorter lifetime $\tau_1$ was fixed at 300ps and the concentration ration $k_1/k_2$ was set as 1. The result (cf. Fig. \ref{fig:fluo_error_flim}) is in accordance with the conclusion from the literature that the optimized gates separation is around 3 times than the shorter lifetime, and when the difference between the two lifetimes increases, a smaller error is achieved \ref{}.
%
%\begin{figure}
%\centering
%\includegraphics{fluo-9}
%\caption{Propagation error on $k1/k2$ for different lifetime pairs and gate seperation.}
%\label{fig:fluo_error_flim}
%\end{figure}
%
%
%\begin{table}[!h]
% \tabcolsep 0pt \caption{optimal gates for pure decay gate selection } \vspace*{-12pt}
%\begin{center}
%\def\temptablewidth{\textwidth}
%{\rule{\temptablewidth}{1pt}}
%\begin{tabular*}{\temptablewidth}{@{\extracolsep{\fill}}ccccccc}
% $N_{gates}$ & $N_{gate combinations}$ & Optimization time & \multicolumn{2}{c}{Optimized gates error} & \multicolumn{2}{c}{Evenly spanned gates error}\\
%%  &$\tau_1$&$\tau_2$&$\tau_1$&$\tau_2$&$\tau_1$&$\tau_2$\\
%
%\hline
%2 &  7140     & &100\%, 18\%     &23.08\% &\\
%3 &  280840    & &100\%, 18\%, 16\%  &17.93\% &\\
%4 &  8214570   & &100\%, 18\%, 17\%, 16\%  &15.40\% &\\
%5 &  190578024  & &  &   &\\
%... &        & &  &   &\\
%118 & 7140     & & &4.38\%  &
% \end{tabular*}
% {\rule{\temptablewidth}{1pt}}
% \end{center}
% \label{opti_para}
% \end{table}
%
%In tomographic settings, the exponential decay is altered by the diffusing media, thus A and B are replaced by the convolution of the excitation tpsf and the pure decay curve in time. The procedure of this method is:
%\begin{itemize}
% \item Obtain excitation time-gated measurements
% \item Convolve the measurements with decay curves for different lifetimes
% \item Calculate the optimum gate sets for a few $k$ ratios that are falling in a biologically reasonable range.
% \item Obtain fluorescence signals at the target gates
% \item Prepare the time-gated Jacobians with lifetime convolution and reconstruct
%\end{itemize}
%The convolved excitation measurements are not the actual fluorescence signal decay, however it provides estimation of the time profiles of the detected signals.
%
%We then applied the approach on the data from the model described in Section \ref{sec:fluo-simulation}, with measurements having 120 300ps gates. In the error propagation function, the maximum gate has been set to be 1500 counts.
%
%Optimal gate sets using diffused fluorescent decays (average for all detectors)
%2 gates: 12+97;86.65\% error
%3 gates; 13+103+104; 67.15\% error
%4 gates; 14+107+108+109, 58.72\% error
%
%\begin{table}[!h]
% \tabcolsep 0pt \caption{reconstruction result for different gate selection. All the values herein are provided in \%.} \vspace*{-12pt}
%\begin{center}
%\def\temptablewidth{0.9\textwidth}
%{\rule{\temptablewidth}{1pt}}
%\begin{tabular*}{\temptablewidth}{@{\extracolsep{\fill}}ccccccc}
%  Gates position& Expected error (\%) & \multicolumn{2}{c}{Quantif. (\%)}&\multicolumn{2}{c}{Crosstalk (\%)}\\
% &$k_1/k_2$&$k_1$&$k_2$&$k_1$&$k_2$\\
%\hline
%100\%, 60\%   &aa.aa &56.20 &29.66 &36.28 &84.73\\
%100\%, 20\%   &aa.aa &55.27 &33.34 &22.63 &45.54\\
%Optimized 2     &aa.aa &   &   &   & \\
%100\%, 60\%, 20\% &aa.aa &61.42 &38.95 &25.47 &51.72\\
%100\%, 40\%, 20\% &aa.aa &56.78 &42.69 &23.38 &40.60\\
%Optimized 3     &aa.aa &   &   &   &
%\end{tabular*}
% {\rule{\temptablewidth}{1pt}}
% \end{center}
% \label{opti_para}
% \end{table}
%
%In this section, we have evaluated a new method for selecting the optimal gates for FMT based on the error propagation theory. The results establish that similar result can be achieved applying the optimum gates compared to using more gates in the inverse problem, with a reduction of the instrumental acquisition time, the calculation time for the double convolution, and the matrix size to invert. This method is therefore valuable in improving both time and memory efficiency in FMT reconstructions.
%
%\section{FRET application}
\section{Discussion}
In this chapter we proposed a new formulation to compute time-resolved fluorescence Jacobians for FMT. We validated this MC formulation for time-gated datatype to simultaneously image two fluorophores using \emph{in silico} and experimental studies. We validated this MC formulation for time-gated datatype to simultaneously image two fluorophores using \emph{in silico} and
experimental studies. The formulation is computationally efficient and offers three distinct advantages. First, once the excitation simulation is computed, it can be used to rapidly calculate the weight functions for multiple fluorophores with different lifetimes. Thus this method becomes attractive for lifetime multiplexed studies. Second, this scheme can be extended to multi-spectral imaging where the absorption coefficients can be modified in (\ref{eq:W_prime_3}) to accommodate the changes in extinction coefficients at different wavelengths without the need of simulating another set of absorption coefficients. Third, using (\ref{eq:W_prime_1}) the difference between $\mu_a^x$ and $\mu_a^m$ can be taken into account, noticing that there are no significant increase of computational burden using (\ref{eq:W_prime_1}) and (\ref{eq:W_prime_3}) thanks to the WMC method~\cite{Alerstam2008a}. Such flexibility in incorporating differences in spectral absorption allows its use for compounds with relatively large Stokes shift (e.g. quantum dots). Here we investigated a 20\% change in the \emph{in silico} study and an improvement less than 5\% was achieved using the corrected $\mu_a^m$, implying the proposed approach is robust for reconstructions in the
NIR window. On the other hand, the scattering coefficient cannot be rescaled like the absorption coefficient due to the transmittance geometry employed in this work~\cite{Alerstam2008a}. While an inaccurate estimation of the scattering coefficient will result in errors in the reconstruction, $\mu_s$ varies more slowly (seldom more than 20\%) over the spectral range in the NIR window~\cite{Mourant1998}. Previous studies have reported that 20\% change in $\mu_s$ will result in less than 15\% error in diffuse optical tomography~\cite{Cheng1999} and similar results are expected in FMT.

We validated the proposed formulation by reconstructing two different compounds simultaneously based on one excitation/emission wavelength pair without using an lifetime based unmixing algorithm on the diffuse measurements. We employed the time-gated datatype to retain the appeal of time-resolved measurements, improved resolution and lifetime contrast unmixing. While using early gates increases the resolution, the crosstalk between fluorophore is significant due
to their equal contribution to the signal when used alone. The introduction of late gates in reconstruction is required to improve separation and quantification of fluorescent probes with lifetime contrast. The combination of early and late gates allows one to trade off the advantages and disadvantages of different time gates and thereby obtain enhanced resolution and crosstalk simultaneously. Employing 5 gates, the fluorescent inclusions were accurately resolved in the region of interest with minimal crosstalk. However, not all the individual metrics are optimal as an emphasis was put on reduction of crosstalk in this work due to the particular application considered. In the case of applications with only one compound, it will be probably possible to identify a set of gates that allow superior resolution and quantification~\cite{Niedre2010}. Moreover, a reconstruction scheme, in which the spatial a priori information derived from an early gate reconstruction is incorporated in reconstruction based only on late gates (or even CW datatype for enhanced SNR), could provide superior results. Also, in the experiment conducted herein, the two compounds were spatially separated. The method proposed here however does allow the reconstruction of the concentration of two fluorophores that are co-localized. This cannot be achieved when employing methods that reconstruct both quantum yield and lifetime. And this provides a new method to monitor the multiple markers \emph{in vivo}.

\section{\label{}Summary}
In this chapter we developed a new formulation for performing time-resolved FMT based on the Monte Carlo forward model. The derived formulation allows to compute fluorescence Jacobians within the fluorescence Born formulation in a computationally efficient manner. Based on standard assumption in the field of FMT, i.e. that the optical properties are similar at the excitation and emission wavelength, the time-resolved fluorescence Jacobians for multiple fluorophores can be computed from one excitation simulation. We investigated the use of this new formulation in the context of simultaneous tomographic imaging of multiple fluorophores without using any unmixing algorithm. We considered the use of time-gated datatypes spanning the full TPSFs, including the early gates for which the diffusion equation is known to be inaccurate when compared to the MC method. The \textit{in silico} and experimental studies demonstrated that this formulation can simultaneously reconstruct two fluorophores with improved resolution and reduced crosstalk when multiple gates are employed, whereas the two fluorophores cannot be resolved if only early gates are employed. Moreover, we demonstrated the successful use of this model combining with structured illumination strategies in the time domain.
